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ABSTRACT 

Magnetohydrodynamic simulations of fully convective, rotating spheres with volume heating near 
the center and cooling at the surface are presented. The dynamo-generated magnetic field saturates 
at equipartition field strength near the surface. In the interior, the field is dominated by small-scale 
structures, but outside the sphere by the global scale. Azimuthal averages of the field reveal a large- 
scale field of smaller amplitude also inside the star. The internal angular velocity shows some tendency 
to be constant along cylinders and is "anti-solar" (fastest at the poles and slowest at the equator). 

Subject headings: stars: low-mass, stars: pre-main sequence, stars: magnetic fields, MHD, convection, 
turbulence 



1. INTRODUCTION 

The Hayashi track in the Hcrtzsprung-Russell diagram 
characterizes young stars in hydrostatic equilibrium that 
are fully convective. Other fully convective stars are 
low-mass main sequence stars (M dwarfs), and some 
cool giants. These stars show strong magnetic activity 
as is evidenced by chromospher ic emission in Ha (e.g. 
!Hawlevlll993t lHawlev et al.lll999]) and by Zeeman broad- 
ening of classical T Tauri stars (e.g. Uohns-Krull et all 
1999b). In the latter case, the stars are generally rapidly 
rotating with rotation periods of just a few days, and 
it is known that the magne tic field shows strong de - 
partures from axisymmetry ij.Tohns-Krull et aLl ll999aL 
However, for less massive stars (M9 dwarfs and beyond) 
there is a sharp decline i n chromospheric magnetic activ- 
ity fe.g. lGizis et alJl2000|L which may be connected with 
dust formation and the a lmost fully neutral photospheres 
(jMohantv fc Basrill2nml) . 

Despi te some progress in lo w-resolution Doppler imag- 
ing (e.g. l.Toncourt et aTlll994l) . not much is known about 
the surface differential rotation of these stars, and even 
less is known about their internal angular and meridional 
velocities. Theory suggests that the absolute differential 
rotation in fully convective stars decreases with increas- 
ing overall angular velocity due to rotational quench- 
ing of the turbule nt transport effect that causes the dif- 
ferent ial rotation (iKjjke relalJ [19_9_3j; iKitchatinov et all 
119941: iKiiker fc Riidigerl 119971 11999ft . As in the solar 
case, the equator is still predicted to rotate more rapidly 
than the poles. However, some observations of rapidly 
rotating stars support what is sometimes referred to 
as "anti-solar" differential rotation, where the equator 
spins less rapidly than the poles llBarnes et aD 120041: 
IKitchatinov fc R.iidigerl 12001 1 Weber et all 120051) Since 
differential rotation enters as an important ingredient in 
dynamo theory, it is important to develop self-consistent 
models of the large-scale velocity field in fully convective 
stars. 

Magnetic field generation in fully convective stars is 
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also interesting from a dynamo theoretical point of view. 
With the realization that the magnetic field inside stars 
might be highly intermittent and concentrated in thin 
flux tubes, the question of storing such intermittent 
and strongly buoyant magnetic fields over the course 
of the 11 year cycle became a growing concern (e.g. 
iMoreno-Insertisl Il983|) . This led to the proposal that 
dynamos in convective shells (as in the case of the Sun) 
might operate at or below the bottom of the convec- 
tion zone. This scenario would not be applicable to fully 
convective stars, because they lack the overshoot layer 
where strong flux tubes could be stored. However, it 
is known that the chromospheric activity does not dis- 
appear for l ater spectral types, i.e. towards fully con- 
vectiv e stars (|Vilhulll984tlVilhu et al 1119891 iBerger et, al l 
2005). It has theref ore been claimed l|Durnev et, alll993t 
lHawlev et alJ2000l) that fully convective stars lack large- 
scale magnetic fields, but can still have small-scale fields 
generated by non-helical near-surface turbulent dynamo 
processes. 

Attempts to model s uch small-scale dynamo action 
ijDorch fc Ludwigl 120021 have however led to the con- 
clusion that the photospheric conductivities of M-dwarfs 
are most probably too low to allow for local small-scale 
dynamo action This would imply that the observed mag- 
netic activity must be due to dynamo action in deeper 
layers. 

From a kinematic m ean-field a 2 dynamo model, 
IKiiker fc Riidigerl l)1999fl predicted that rapidly rotating 
(Coriolis number of 3 or larger) fully convective stars 
generate a non-axisymmetric steady magnetic field of 
quadrupolar symmetry and azimuthal order m = 1 that 
looks roughly like a dipolc field with the dipole axis lying 
in the equatorial plane. 

Global models of convective dynamos are still in 
their infancy, even though some t remendous pro gress 
was made some 2 years ago when iGilmanl l)1983|) and 
iGlatzmaierl (|1985|) presented the first simulations of dy- 
namos in a spherical shell representing a solar-like con- 
vection zone. These models predicted cyclic magnetic 
fields propagating toward the poles, in contrast to the 
solar case. The reason for this discrepancy remains a 
matter of debate even today, when much higher numer- 
ical resolution is available. Recent simulations still pre- 
dict angular velocity to be roughly constant on cylinders, 
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although some simulations show at least a tendency to- 
ward solar-like angular velo city contours ( Micsc h et alJ 
120001 iBrun fe Toomrel l2002). Recent simulations of dy- 
namo action in spherical shells now begin to produce 
useful models of global turbulent dynamos l)Brunll2004l 
IBrun et al.ll2004j) . Meanwhile, such glo bal models have 
also been applied to core convection (Brownine:^^aU 
2QQ4) and to dynamo action in these cases IjBrun et all 



2005) 



In this paper, we present global dynamo simulations 
in spheres using a Cartesian grid, i.e. the sphere is em- 
bedded in a cubic box. This may seem to be an unnat- 
ural approach to spherical geometry, but it has distinct 
practical advantages. First, it avoids the coordinate sin- 
gularity at the center when using spherical coordinates, 
without invoking expensive transformations from and to 
spherical harmonics. Second, this approach has proven 
useful in view of computational simplicity and numerical 
parallelization efficiency; it has recently been applied by 
a number of groups to purely hydrodyn amic simulations 
UPorter et al 1200(1 iFrevtag et~atll20?rl I Wood ward et a.U 
2003), and attempts have alrea dy been mad e to model 
dynamo action in this approach l|Dorchll2004[) . 

2. THE MODEL 

2.1. Basic setup 

In our model the star is described as a spherical subre- 
gion of radius R of a cubic box of size L^ ox . The gas in 
the box is governed by the usual equations of magnetohy- 
drodynamics (see below) with impenetrable boundaries 
on the box faces such that the mass M^ ox in the box is 
conserved. If the gravitational well <i>(r) is sufficiently 
deep, most of the mass M of the star is concentrated 
near the center, soM« Mb ox - Using a Newtonian cool- 
ing term in the energy equation, the temperature outside 
the star is kept close to the nominal surface temperature 
of the star, T sur f . An entropy gradient is maintained by 
prescribing a distributed energy source TL(r) at the cen- 
ter (here r is the spherical radius). The total luminosity 

is then given by L — A-k J H.(r)r 2 dr, and corresponds 
to the energy produced by nuclear burning. We recall 
however, that some young stars on the Hayashi track 
have not ignited yet, and are sustaining their energy 
losses by contraction, which results in a less localized 
energy source than nuclear fusion reactions. Although 
the mass distribution can change during the evolution of 
our model, we have chosen to ignore self-gravity. 

The model is governed by five main input parame- 
ters: mass M, radius R, luminosity L, surface tempera- 
ture T sur f, and average angular velocity f^o- We choose 
parameters that are typical of M dwarfs, but we limit 
the degree of stratification to values that are numeri- 
cally more feasible by choosing a surface temperature 
that is much higher than for real M dwarfs. We also 
keep the Kelvin-Helmholtz time scale at a much smaller 
multiple of the dynamical time scale than what is realis- 
tic. As is common in deep convection simul ations (e.g., 
iChan fc Sofiall98?llBrandenburg et a l. 2005), we do this 
by choosing a luminosity that is much larger than the 
stellar value, and at the same time we keep the radiative 
diffusivity much larger than in reality. Since the Rayleigh 
number is, for a given Prandtl number, inversely propor- 
tional to the square of the radiative diffusivity, a large 



luminosity translates to a small Rayleigh number. The 
restriction to moderate values of the Rayleigh number 
is a common problem of all astrophysically meaningful 
convection simulations. 

Our initial state is derived from a spherically symmet- 
ric, isentropic reference model; for details see Appendix 
lAl This state is perturbed by adding weak velocity and 
magnetic fields that are both random. 

2.2. Equations 

In the computational domain, — Lbox/2 < x,y,z < 
-^box/2, we solve the equations of compressible magneto- 
hydrodynamics , 

Ding 



Di 
Du 
"57 



Vu, 
Vp j x B 



g g 
- V$ - 2fl x u + f d 

— = u x B - r/^oJ , 



- V • (qvS) 
Q 



Mo?7J 



2gvS 2 , 



(1) 

(2) 
(3) 

(4) 



where g and p denote mass density and pressure of the 
fluid, s and T are specific entropy and temperature, u is 
the fluid velocity, v the kinematic viscosity, $ the gravity 
potential, fio the angular velocity of the reference frame, 
fd is an artificial damping force discussed in Sec. 12.31 and 
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is the traceless rate-of-strain tensor. The magnetic vector 
potential A is related to the flux density B = V x A and 
the current density j = V x B//io, and ?7 denotes the 
magnetic diffusivity. Volume heating H. and cooling C are 
described in Sec. 12.31 below. The radiative conductivity 
K is related to the thermal diffusivity x = K/(c p g). In 
the numerical calculations shown below, we assume \i v ^ 
and rj to be constant across the whole box. Our equation 
of state is that of a perfect gas with adiabatic index 7 = 
5/3. 

For the gravity potential <&(r) we choose a Pade 
approximation obtained from our isentropic reference 
model (see Appendix |A"|) . 
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with r' = r/R; we find that the coefficients, ao = 2.34, 
a 2 = 0.44, a 3 = 2.60, b 2 = 1.60, b 3 = 0.21, yield a good 
approximation both inside and outside the star. 

Note that, while retaining the Coriolis force term, we 
neglect the centrifugal force. This is necessary for prac- 
tical reasons, since together with the luminosity our tur- 
bulent velocities w rms are exaggerated and we thus need 
far too large angular velocities in order to reach realistic 
Coriolis numbers [see Eq. Ijl3|l below]. It would there- 
fore be unrealistic to include the strongly exaggerated 
centrifugal force in the expression for the inertial forces. 
We emphasize that this kind of restricted mechanics does 
not violate the balance of angular momentum L in any 
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Fig. 1. — Comparison of luminosity function L r (r) = 
Jg H(r') 47rr' 2 dr' according to Eq. IA4I (dashed line) with our 
Gaussian parameterization 0. The choice u>l = 0.162 R gives an 
excellent fit, while the narrower profile would be more appropriate 
for a heavier star. 

significant manner: The component L z parallel to the 
rotation axis is strictly conserved (the centrifugal force 
is a central force for that axis), while the other two com- 
ponents are small for a nearly axisymmetric system. 

Our boundary conditions on the faces of the cubic box 
are impenetrable free-slip conditions for velocity (u n = 0, 
c? n Utan = 0), and normal-field conditions for the mag- 
netic field {d n B n = 0, Bt an = 0). 

All numerical calculations were done using the Pen- 
cil Code, 4 a high-order centered finite-difference code 
(sixth order in space and third order in time) for solving 
the compressible hydromagnetic equations. Weak shock- 
capturing viscosities were used to cope with localized, 
transient events of supersonic flow. A high-order upwind 
scheme is used for the advection operators for density 
and entropy; see Appendix FBI 

2.3. Profile functions 

As outlined in Sec. 12.11 the thermal structure of the 
star is maintained by prescribing a certain distribution of 
heating and cooling functions inside and outside the star, 
respectively. The profile functions depend on spherical 
radius r = {x 2 +y 2 +z 2 ) 1 / 2 . In the exterior, r > R, we 
add a velocity damping term in order to prevent excessive 
velocities outside the star, which are not directly relevant 
to the dynamics inside the star. 

The central parts of the sphere are heated according 
to a normalized Gaussian profile, 



H(r) = 



L 



(7) 



(2irw 2 L ) 3 / 2 

which gives an ex celle nt fit to the heating rate calculated 
according to Eq. l|A4jl for our isentropic reference model 
if the width wl of the nuclear burning region is chosen as 
wl = 0.162 R. Fig.^shows a comparison of the resulting 
luminosity from the two parameterizations. Most of our 
simulations use that value of wl, while some runs have 
been carried out with u>l = 0.1 R, which would be more 
appropriate for a more massive star. 



4 http://www.nordita.dk/software/pencil-code — This code 
uses the Message .Passing interlace (MF1) library for communi- 
cation between processors and runs quite efficiently on clusters. 
Toroidal averages, spectra, and other diagnostics can be calculated 
during the run, which avoids extensive post-processing of the data. 



For r > R, we apply a Newtonian cooling term of the 
form 

- C{r) = -gcp — fext(r) (8) 

7"cool 

to keep the temperature close to the surface value T sur f . 
Here, / ex t(?") is a profile function that smoothly interpo- 
lates between for r <C R and 1 for r ^> R. Our profile 
function is a tanh profile, 



/extO") = \ { 1 + tanh ^ c ° Ql 

2 V Wcool 



(9) 



where i? C ooi and w coo \ denote the position and width of 
the transition. We have chosen w coo \ = 0.05 R, and 
-Rcooi = 1.05 i?, i.e. slightly larger than the stellar ra- 
dius, in order to reduce the influence of the cooling term 
ijHJ inside the star. In the present model, the exterior 
has practically constant temperature (= T sur f), i.e. no 
attempt is made to model the hot corona of the star. 
In fact, since we have to restrict ourselves to moderate 
stratification, the temperature ratio between the center 
and the surface of the model is less than 10. 
Outside the star, a damping term 



fd = fext(r) 



(10) 



is applied in the equation of motion to limit flow speeds 
to moderate values while still allowing the exterior to 
react to sudden disturbances from the stellar surface with 
sufficient flexibility; the profile / xt(f) is the same as for 
the cooling term, i.e. Eq.© with = 0.05 R, and — 
1.05 R. By imposing fixed radial profile functions for 
surface cooling and velocity damping, we suppress the 
possibility of i rregular surfaces tha t would develop, e.g. 
in red giants l(Frevtag et alJl2002T) . but this would not 
apply to real M dwarfs. 

2.4. Dimensionless parameters 

As mentioned in the beginning, our model is governed 
by the five basic input parameters: M, R, L, T SUI f, and 
CIq, From these, we can construct three dimensionless 
quantities that characterize our model: the stratification 
parameter 



rfh 



GM/R 



(11) 



(where c SjSur f is the sound speed at the surface and G is 
Newton's gravity constant), the dimensionless luminosity 



C = 



L 



y/G 3 M 5 /R 5 ' 



(12) 



and the Coriolis number (or inverse Rossby number) 



Co = 2fl Q R/u, r 



(13) 



where u rms is the root-mean-square velocity based on 
a volume average over the full sphere. The remain- 
ing degrees of freedom determine the natural units of 
our system. In particular, length will be measured in 
units of the stellar radius [x] — R, velocity in units 

of [u] = y/GM/R, density in units of [g] = M/R 3 , 
and specific entropy in units of [s] = c p . This implies 
that time is measured in units of the dynamical time 
[t] = {GM/R 3 )- 1 / 2 and the magnetic field is measured 
in units of [B] = ^fj, [e] M = VJM)G M/R 2 . 
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TABLE 1 

Summary of runs discussed in the paper. 



Run 


rcsol. 


V 












c 


n 


X 


^kin / sat 


osat 
rms 


Arms 
"rms 


Re 


Rm 


la 


128 3 


6 x 10" 


-4 


4 x 10- 


-4 


3 x 10- 


-4 


0.02 


0.2 


0.017 


0.173/0.164 


0.020 


0.12 


273 


547 


lb 


256 3 


4 x 10" 


-4 


3 x 10- 


-4 


2 x 10" 


-4 


0.02 


0.2 


0.043 


0.184/ 


0.028 


0.18 


388 


775 


lc 


256 3 


4 x 10" 


-4 


3 x 10- 


-4 


2 x 10" 


-4 


0.01 


0.2 




/ 0.130 


0.023 


0.18 


325 


650 


2a 


128 3 


8 x 10" 


-4 


8 x 10" 


-4 


4 x 10- 


-4 


0.02 


0.0 


0.009 


0.239/0.233 


0.018 


0.08 


291 


583 


2b 


128 3 


8 x 10" 


-4 


8 x 10" 


-4 


4 x 10" 


-4 


0.02 


0.5 


0.017 


0.213/0.185 


0.046 


0.25 


231 


163 


2c 


128 3 


8 x 10" 


-4 


8 x 10" 


-4 


4 x 10- 


-4 


0.02 


2.0 


0.021 


0.158/0.129 


0.068 


0.53 


161 


323 


2d 


128 3 


8 x 10- 


-4 


8 x 10- 


-4 


4 x 10- 


-4 


0.02 


5.0 


0.036 


0.112/0.087 


0.099 


1.14 


109 


218 


2c 


128 3 


8 x 10- 


-4 


8 x 10- 


-4 


4 x 10- 


-4 


0.02 


10.0 


0.038 


0.086/0.068 


0.105 


1.54 


85 


170 



Note. — Diagnostic quantities listed arc: kinematic growth rate A of the magnetic field; root-mean-square values 
of velocity and magnetic flux density u rms , 6 rmB ; ratio B Tms /u Tuls for the saturated state; kinetic Reynolds number Re 
(based on it rms ); and magnetic Reynolds number Rm (based on u rms ). For all runs shown here, the star is embedded in 
a cubic box of size Lb ox — 3 R. The gaps for Runs lb and lc are due to the fact that we have not extended Run lb into 
the final saturated regime, but rather lowered the value of C and continued it as Run lc. 



Note that £ is the ratio of the pressure scale height 
at the stellar surface to the stellar radius, so £ controls 
the amount of stratification. The second dimcnsionless 
parameter, C, is the ratio of acoustic (or free-fall, or dy- 
namic) time scale to the Kelvin-Hclmholtz time. For 
realistic models, both £ and C are much less than unity. 
Using typical values for an M5 dwarf (M — 0.21 M©, 
R = 0.27 R Q , L = 0.008 L Q , and T surf = 4000 K), we 
find £ = 2.2 x 10~ 4 and C = 2.4 x 10~ 14 . In the sim- 
ulations presented below, we are only able to reach val- 
ues of £ and C that are somewhat below unity. In all 
models presented here, we have £ = 0.19; for most mod- 
els we choose C = 0.02 (i.e. ~ 10 12 times higher than 
for a real M5 dwarf), while we have C — 0.01 in one 
of the higher resolution runs. The necessity of exagger- 
ated luminosities in nu merical simulat i ons of convection 
was first pointed out bv lChan fc Sofia! (^86). For lower 
values of C, yet higher numerical resolution would be 
required to get sufficiently vigorous convection and dy- 
namo action. 

Other important dimensionless parameters are the 
kinematic and magnetic Reynolds numbers, 

Re = — and Rm = — , (14) 

V 7] 

where U is the root-mean-square (rms) velocity within 
the sphere of radius R. In the present simulations, Re 
and Rm are in the range 100-780 (see Table^). Realistic 
values of the fluid and magnetic Reynolds numbers are 
much larger than what can be achieved in this type of 
simulation. 

3. RESULTS 

The parameters for the runs discussed and presented 
in this paper are summarized in Table ^ Throughout 
this paper, overbars denote azimuthal averages. The rms 
values listed here are also averaged in time. In Runs la- 
ic, luminosity and resolution have been varied, while in 
Runs 2a-2e we have varied the angular velocity ^o- 

The simulations were typically run for about 8rohm> 
where rohm = R 2 /^ 2 ?]) is the diffusive time scale for 
a structure of wave length 2 R. One exception was the 
higher-resolution runs lb and lc, which were only run 
for about 0.3 Tohm each. In all cases, the saturated state 



of the magnetic field was well established (with the ex- 
ception of Run lb, which we did not run long enough) 
and quasi-stationary behavior was reached. To ensure 
that we are not missing any slow trends, we continued 
Run 2c until 12rohm) but found nothing new during this 
somewhat prolonged saturated calculation. 

For all runs listed in Table ^ the box size was Lbox = 
3 R. To investigate the role of the boundaries of the 
numerical box, we did a reference run in a larger box 
(^box = 5 R) at comparable resolution. The results were 
fully compatible with L D ox = 3 R. 

3.1. Radial stratification 

Figure [3 shows density, squared sound speed (propor- 
tional to temperature), Mach number and specific en- 
tropy as function of radius for Run lc. Density and 
squared sound speed vary by a factor of about 5 from 
the center to the stellar surface. Apart from a few lo- 
calized transients, the maximum Mach number is below 
unity and there is no evidence for shocks. The total vari- 
ation in specific entropy is about 0.6 c p . Even in the bulk 
of the convection zone (0.15 < r/R < 0.85) the specific 
entropy has a standard deviation of about 0.05 c p , which 
is still much larger than what mixing-length theory pre- 
dicts for this type of star. This is related to the high 
value of C that we are using, which is also the reason for 
the enhanced entropy values in the core. The location 
of the specific entropy minimum is at r/R 0.93, i.e. 
somewhat below the nominal surface of the star. This is 
primarily a consequence of the rather large width of the 
profile functions for cooling and velocity damping, which 
affect the interior already inside r — R. At that effective 
radius, we naturally get a thin overshoot layer (as found 
in real stellar chromospheres). 

3.2. Hydrodynamic flow patterns 

Since the initial magnetic field is weak (several orders 
of magnitude below saturation), the kinematic phase of 
the dynamo represents the hydrodynamic flow pattern 
in a nonmagnetic scenario. Figure shows an equato- 
rial section of entropy and density for Run lb. One can 
clearly distinguish narrow cool structures (downdrafts) 
that are familiar from box simu l ations of compressible 
convection (e.g. iHurlburt et al.l Il986t iNordlund et alJ 





-1.0 -0.5 0.0 0.5 1.0 

xIR 



Fig. 3. — Equatorial section of entropy (color coded: dark for low, 
bright for high entropy) and density (isolines) during the kinematic 
stage of Run lb (at t = 300 [t] ) , where the magnetic field does not 
affect the dynamics. 



1992) The flows are far from being laminar, as can also be 
seen in the inset to Fig.^J However, given the numerical 
resolution, only a limited range of scales can be resolved, 
as can be seen from the magnetic energy spectrum during 
the kinematic dynamo phase (see next section). 
Figurc[S]shows at-cp average of kinetic helicity u • V x u 



for the kinematic dynamo phase. As expected from the 
action of the Coriolis force on expanding upflows and 
contracting downflows, the helicity is predominantly neg- 
ative in the northern and positive in the southern hemi- 
sphere. If kinetic helicity is connected to a turbulent 
electromotive force, we find a distribution of the a effect 
that is reminiscent of classical mean-field dynamo mod- 
els (e.g. lRobertsll972|) . It should hence not be surprising 
if the flow generates a large-scale magnetic field. 

3.3. Dynamo action 

The turbulent kinetic energy quickly reaches a sta- 
tistically steady state after about 5 dynamical times 
(t w 5 [t] ), while the energy of the initially random mag- 
netic field decays at first; see Fig. This is because 
most of the magnetic energy of the random field is in the 
small scales and thus gets quickly dissipated. The mag- 
netic field then grows exponentially with a growth rate 
A = dlnB rms /di of about 0.04/[t] (for Run lb). 

During the kinematic stage of the dynamo, the mag- 
netic field grows exponentially with the same rate at all 
wavenumbers, so the spectrum remains shape- invariant, 
as can be seen in Fig. [5] The maximum of the magnetic 
spectrum is around k as 3 x 2n/R. The growth time 1/A 
is about one order of magnitude shorter than the global 
diffusive time scale rohm, which is a manifestation of 
turbulent magnetic diffusion. 

At later times, magnetic energy saturates first at the 
smallest scales, while the large scales still accumulate 
energy. Eventually all scales are saturated, but now the 



6 



Dobler et al. 



1 1 1 1 1 1 


M rms 






J 




_ 


J 






/ j1.0 

/ E 


















r 


20 40 60 80 


i i i ■ ■ ■ 


f/W 





100 200 300 400 500 600 700 

t/[t] 



Fig. 4. — Evolution of root-mean-square velocity u TInB and mag- 
netic field B rms (represented as an Alfven speed using g = 0.4 [g] 
to make the two curves comparable) for Run lb. Time is measured 
in units [t] and velocities in units [u] (see Sec, 12, 4^ , The inset shows 
the maximum velocity M max (i) during the onset of convection — 
note the irregular time-behavior. 
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Fig. 5. — Average of kinetic helicity for Run 2c during the 
kinematic phase. Shown is the azimuthal average of u ■ V X u, 
averaged in time from t = 100-300 [t]. 
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Fig. 6. — Spectra of the magnetic field for times t = 
(300, 600, 900, 1200, 1500, 2100) [t] of Run 2c (for this run, expo- 
nential field growth levels off around t 700 [t] ) . Magnetic energy 
increases with time and eventually reaches saturation. At late 
times the largest scales dominate. 



magnetic spectrum peaks at a larger scale than during 
the kinematic stage. As the magnetic field reaches sat- 



io- 3 F 




kRIln 

Fig. 7. — Magnetic (solid line) and kinetic (dashed line) power 
spectra for the saturated phase of Run 2c. While the velocity spec- 
trum peaks around k 2n/R, magnetic energy is distributed more 
flatly around the largest scales. For comparison, the dotted line 
shows the kinetic spectrum during the kinematic dynamo phase. 



uration, the kinetic energy of the flow is decreased by a 
certain amount that depends on the relative importance 
of rotation (see TablcQ. For slowly rotating spheres the 
kinetic energy decreases by only about 10-20% (Runs la 
and 2b), but for more rapidly rotating spheres, where 
the magnetic energy is also much larger, the suppres- 
sion of the kinetic energy is about 50-60% (Runs 2c-2e) . 
The strong dependence of the kinetic energy on the mag- 
netic field strength suggests that the flows are probably 
not strongly turbulent and still governed by a large scale 
more laminar flow pattern. 

Increasing the resolution by a factor of 2, while at the 
same time decreasing dissipative effects (cf. Runs la, lb), 
we see that the growth rate increases significantly (by a 
factor of 2.5, see Tabled, but in the saturated state the 
rms velocity changes insignificantly. The rms magnetic 
field increases by about 40%, which is rather large and 
may be a consequence of the dynamo not being strongly 
supercritical. Decreasing the luminosity by a factor of 
2 (cf. Runs lb, lc) decreases rms velocity and magnetic 
field only by about 20%. 

Figure shows spatial spectra of kinetic and mag- 
netic energy. Kinetic energy peaks at a wavenum- 
ber of about k p w 1 x 2ir/R 61a:]" 1 , which cor- 
responds to the energy-carrying scale [smaller scales 
have a negligible contribution to the total kinetic en- 
ergy J E(k)Ak\. The corresponding turnover time is 
t = (frms^p) 1 ~ 2[i]. This results in a normalized 
growth rate At = 0.08, which is comparable to the 
values for both helically and nonhelically forced turbu- 
lence simul ations where At = 0.0 3-0.1, see iBrandenburel 
({20011) and lHaugen et alJ l|2004D . respectively. The sat- 
uration value of magnetic energy is typically an order of 
magnitude below the kinetic energy of the turbulence 
for the slowly rotating models. This is quite similar 
to the ratio found in earlier simulations of convection 
driven dynamos in Cartesian and spherical geometries 
iMeneguzzi fc Pouauetll989HNordlund et alTl 992: Brun 
2004). For the faster rotating Runs 2c-2e, however, mag- 
netic energy is comparable to kinetic energy. 

Theoretically, there is always the possibility of differ- 
ent solutions to a nonlinear problem, depending on the 
initial conditions. This possibility has been anticipated 
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in co nnection with the geodynamo (Ro berts fc Sowardl 
1992), but it has so far not been seen in any turbulen t dy- 
namo simulation (e.g. iGlatzmaier fc Roberts! ITool . In 
principle there is even the possibility of so-called self- 
killing dynamos that decay after full saturation has been 
reached, but such behavior has so far only been found 
under rathe r artificial conditi ons and in the absence of 
turbulence l)Fuchs et alJ 11999^1 . In some cases we have 
restarted our simulations from a snapshot that has been 
obtained for different parameters. In such cases we al- 
ways recovered statistically the same solution that was 
obtained in the standard way by starting from a weak 
seed magnetic field and a non-convecting initial state. 
This was further confirmed by restarting Run 2c with 
a 10 10 times weaker initial field, which led to an indis- 
tinguishable time history of u max for t < 5 [£], 5 and to 
statistically equivalent behavior for larger times. 

The saturation appears to happen on a dynamical time 
scale, i.e. we see no evidence for resistively limited sat- 
uration, as it was found in h elically forced simula tions 
in a triply periodic domain l|Brandenburd l200lh . In 
Run lb, the total simulated time t ss 700 [t] corresponds 
to about 2"rohm- The nonresistive saturation behavior 
could be due to the fact that in the present simulations 
the boundaries are open and permit a magnetic helic- 
ity flux across the equatorial plane and out of th e box 
l|Brandenburg fc Doblerll2001i lBrandenburd l2005). An- 
other possibility is that the magnetic Reynolds number is 
still too small for magnetic helicity conservation to have 
an effect. 

3.4. Dependence on rotation rate 

In Runs 2a-2e, we vary the rotation rate S7o , while 
keeping all other parameters fixed. As the rotation rate 
is increased, the root-mean-square velocity of the tur- 
bulence decreases. This is to be expected, because the 
presence of r otation is known to d elay the onset of con- 
vection (e.g. lChandrasekharlfl961l) . The rms magnetic 
field strength increases monotonically with the Coriolis 
number (or the rotation rate flo) for the runs shown. 

In the absence of rotation (Run 2a), there is no net 
helicity or net shear and hence no reason for the genera- 
tion of a large-scale magnetic field. However, we still find 
that the magnetic energy increases, albeit more slowly 
and to a lower value than for the rotating runs. This is 
a man ifestation of the 'fluctuating' or 'small-scale' dy- 
namo l)Kazantsevlll968HMeneguzzi et al.1ll981tlCattaneol 
1999), which requires a considerably larger value of the 
magnetic Reynolds number than the helical dynamo, and 
there are indications that Run 2a is only mildly super- 
critical. 

Comparing the magnetic energy spectra for runs with 
different rotation rates (Fig. |SJ, we find that the mag- 
netic energy at the large scales increases with Qq at 
least up to Oo ~ 5 (corresponding to Co » 100), while 
the small scales are only weakly affected by rotation. 
The saturation for rapid ro tation has been predicted by 
mean-field dynamo theo ry, ( Rii diger fc Kichatinovlll99'3t 
lOssendriiver et al"ll2001j) . and for even larger Hq, one ex- 
pects a reduction of large-scale dynamo efficiency. How- 
ever, for our models the peak dynamo efficiency occurs 

5 The resulting Lyapunov time scale is a few turnover times 
fi/ti Ims , as one would expect. 
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Fig. 8. — Magnetic energy spectra for the saturated states of 
Runs 2a, 2b, 2c, 2d, and 2e. 




Fig. 9. — Three-dimensional visualization of the magnetic field 
for Run 2c at t = 2600 [t] (saturated phase). Magnetic field lines 
are shown, together with the surface of the sphere. 



for rather large values of Co. Another surprise is that 
dynamo activity at small scales (kR/2n > 10) is not 
quenched for 'superfast' rotation, although the Coriolis 
force should play a significant role until k reaches signif- 
icantly larger values. 

3.5. Large-scale field structure 

Although a lot of the magnetic energy is due to the 
small-scale structures, as can be seen in the magnetic 
energy spectra, outside the star the field shows features 
of a dipole-like structure with a noticeable contribution 
from the first few multipoles; see Fig. [5] where we show 
a visualization of the three-dimensional magnetic field 
lines. 

A more quantitative presentation of the large-scale 
magnetic and velocity fields is obtained by considering 
azimuthal averages, as shown in Fig. llOl for one snapshot 
of Run 2c. The magnetic field shows a clear large-scale 
component with predominantly quadrupolar symmetry 
with respect to the midplane, but still including dipolar 
contributions. The velocity field shows little large-scale 
structure and varies strongly in time. 

We find a considerably more regular structure when 
applying time averaging to the azimuthally averaged 
data. In Fig. ^] we show, for the saturated states, the 
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Fig. 10. — Azimuthal averages of magnetic field (a) and velocity 
(b) for Run 2c at time t = 2100 [t] (saturated phase). Panel (a) 
shows poloidal field lines of the ^-averaged magnetic field (orien- 
tation of the field lines in the top half is predominantly counter- 
clockwise), superimposed on a color/gray-scale representation of 
the azimuthal mean field, B v , with bright colors representing 
Btp > and dark colors representing B v < 0. Note the mixed 
parity of the structure of the mean field with a strong quadrupolar 
contribution. Panel (b) shows vectors of the mean poloidal velocity 
superimposed on a color /gray scale representation of the mean an- 
gular velocity <5f!(r, z) = u^/(r sin 9) with bright colors for <5f! > 
and dark colors for <5f! < 0. Note that for different snapshots in 
time the velocity field u looks very different. 



correspondingly averaged magnetic fields (top row) and 
velocity fields (bottom row) for four different runs with 
rotation rate f^o increasing from left to right. We note 
that for all runs the averaged velocity field changes very 
little from the kinematic to the saturated stage of the 
dynamo. 

With this averaging, we find almost perfect quadrupo- 
lar symmetry for B in Runs 2b and (slightly less pro- 
nounced) 2c. On the other hand, Runs 2d and 2e show 
very pronounced hemispheric asymmetry that appears 
to be relatively long-lived. For Runs 2b and 2c, the ve- 
locity field shows a meridional circulation pattern that 
is directed outwards at the equator, and surfaces of con- 
stant angular velocity u> are approximately cylindrical. 
For the rapidly spinning Runs 2d and 2e, the asymme- 
try in the magnetic structure is reflected in differential 
rotation and meridional circulation. 

An explicit measure for the efficiency of large-scale field 
generation is the ratio 



<7B 



«B) t 2 > 



1/2 



(15) 



which is given in Table [3 A similar quantity q u is 
also defined for the velocity. Here the overbar de- 
notes azimuthal averaging, (-) t represents time averag- 
ing, while (-) rz denotes spatial averaging over the sphere, 

and B ims = ((B 2 ) t ) rz . This ratio is q& — 0.19 for 
Run 2b and increases further with the rotation rate. 
These values are quite large, suggesting that large-scale 
field generation is quite efficient. However, in forced 
turbul ence simulations with open boundaries and no 
shear (|Brandenburg fc Doblerl 12001^ . q& decreases with 
increasing magnetic Reynolds number. On the other 
hand, simulations of forced turbulence suggest that the 



presence of shear is critical for allowing the dynamo am- 
plitu de to be independen t of the magnetic Reynolds num- 
ber (|Brandenburd l2005) . Further numerical simulations 
are necessary to see whether the same behavior occurs 
here as well. 
The ratios 



Pu 



Pb 



_ «^>;> 



2 \ 1/2 



2 \ 1/2 



(16) 



quantify the importance of the azimuthal components in 
the ip- and t-averaged fields. In the non-rotating case, 
p u is very small, indicating that systematic azimuthal 
flows are weak. With increasing angular velocity, how- 
ever, p u at first increases as systematic differential ro- 
tation evolves. The slight decline of p u in Runs 2d and 
2e is connected with the stronger magnetic field in those 
cases. In fact, both the azimuthal component ((u^) t 2 )r/ 2 

and the poloidal component ((m po i) 2 )r/ 2 decrease mono- 
tonically from Runs 2b to 2e because of the increasing 
large-scale magnetic field. 

The ratio pb is 1.1 in the nonrotating case, it has a 
maximum already for Qq = 0.5 (Run 2b) and then de- 
clines. The total magnetic energy continues to grow with 
higher rotation rates (see Table indicating that with 
increasing rotation rate the poloidal field continues to 
grow, while the toroidal field remains roughly unchanged. 
This seems to be in qualitativ e agreement with the sim- 
ulations bv lBrun et alJ l)2004l) for dynamos in convective 
shells without tachocline. On the other hand, this ratio 
is in all cases small compared to what is expected for 
the Sun, where the tachocline can be expected to have a 
strong effect. 

With the exception of Run 2b, the absolute ampli- 
tude Alo of the differential rotation for kinematic and 
saturat ed states is similar. This is co nsistent both with 
theory l)Kitchatinov fc R.iidigerlll999D and with observa- 
tions showing only a weak dependence of surface dif- 
ferential rotation on stellar rotation for late type stars 
l)Barnes et alJ 120051) . which are however not fully con- 
vective. 

Next, we show the energy ratios 



_ «u 2 )*), z 



0"B 



d2 
rms 



(17) 



quantifying the fraction of energy contained in the ax- 
isymmetric part of u and B. For most of the rotating 
runs, (7b increases drastically from the kinematic stage 
to saturation. This is another manifestation of the trend 
towards large-scale fields once the small scales are sat- 
urated (see Fig. EJ. On the other hand, cr u , which is 
strongly reduced by rotation, is not severely affected by 
the magnetic field saturation. 

Finally, we consider the parity of the mean field with 
respect to the equatorial plan e. Earlier work on mean- 
field dynamos in full spheres ijBrandenburg et"aT]ll989j) 
has shown that for weakly supercritical dynamos the par- 
ity is dipolar (antisymmetry with respect to the equator). 
However, as the dynamo becomes more supercritical, the 
parity can become quadrupolar (symmetric with respect 
to the eq uator), but mixed and chaotic behaviors are also 
possible l|Covas et al.lll999|) . Parity of the mean field can 
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Fig. 11. — Azimuthal averages for Runs 2b-2e as in Fig. llOl but now the fields are additionally averaged over 1000—2000 time units 
[t] (about 150-300 turbulent turnover times, or 4-8 rohm) during the saturated state. Angular velocity Qo increases from left to right; 
the panels are labeled by the names of the individual runs. Top row: magnetic field. Orientation of the field lines is predominantly 
counterclockwise in the top hemisphere and clockwise in the bottom hemisphere in all four panels. Bottom row: velocity. Note that the 
angular velocity shows some tendency to be constant along cylinders for Runs 2b and 2c, while magnetic and velocity field are asymmetric 
for Runs 2d and 2e. The amplitudes of magnetic and velocity fields have been scaled individually for each image; absolute values of B rms 
and n rms are given in Table HI 



TABLE 2 

Additional diagnostic quantities for Runs 2a-2e. 



Run 


O 


<3u 


9B 




Pb 


Aoj kin 1 sat 


kin / sat 


kin / sat 
CT B 


P 


2a 


0.0 


0.30 


0.02 


0.04 


1.10 


0.53/0.24 


0.240/0.214 


0.071/0.079 


-0.12 


2b 


0.5 


0.38 


0.19 


1.44 


2.95 


1.23/0.97 


0.222/0.221 


0.077/0.147 


+0.99 


2c 


2.0 


0.14 


0.39 


1.80 


2.36 


0.56/0.22 


0.062/0.079 


0.124/0.271 


+0.91 


2d 


5.0 


0.14 


0.40 


1.59 


1.32 


0.82/0.26 


0.054/0.064 


0.046/0.271 


+0.28 


2e 


10.0 


0.14 


0.46 


1.27 


1.17 


0.50/ 0.14 


0.080/0.063 


0.047/0.326 


-0.28 



Note. 



All quantities refer to the saturated state, unless explicitly labeled with 'kin'. 



be quantified as 



E S -E A 
E S + E A 



(18) 



where Es and E\ denote the energies contained in the 
symmetric and antisymmetric parts of B. The values 
of P are listed in Table El for the saturated stage of the 
dynamo. For Runs 2b and 2c the mean field is nearly 



quadrupolar P w +1, as is also evident from Fig. II II 6 
Both in the absence of rotation and for strong rotation 
the mean fields are of more mixed parity character. Com- 
parison with mean-field dynamos would suggest that our 

8 Note that the mean field for the snapshot shown in Fig. [5] 
is still predominantly symmetric, although its poloidal field lines 
look quite dipolar. This is due to three factors: the dominance of 
the azimuthal component, the non-axisymmetry of the field, and 
a hemispheric asymmetry that is not very prominent in Fig. |3 
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Runs 2b and 2c are in the "more supercritical" regime. 
However, we have not found a case that would clearly 
belong to the weakly supercritical regime, where dipo- 
lar fields are expected. In addition, mean-field theory 
would suggest cyclic mean fields, that have also not been 
seen. It is possible that such features would emerge in 
a direct simulation only after many more turnover times 
than what has been possible here. 

The ratios q and p allow us primarily to assess the 
mode of operation of the dynamo found in the simu- 
lations. In particular, they provide a sensitive tool to 
assess the possible dependence of the magnitude of large- 
scale fields on the magnetic Reynolds number. Observa- 
tionally, of course, only the limit of very large magnetic 
Reynolds numbers is relevant. More detailed comparison 
of the q and p ratios with observations is hampered by the 
fact that the magnetic field in the star's interior cannot 
be measured with current techniques. The interpreta- 
tion of proxies such as filling factors and the appearance 
of magnetic fields at the stellar surface may be prema- 
ture as long as we do not fully understand the connection 
between magnetic fields at the surface and the interior. 
For example, the interpretation of bipolar spots in terms 
of distinct flux tubes may not be valid and hence the 
presence of spots may not necessaril y indicate a high de - 
gree of intermittency in the interior l)Brandenburel 2005 ) . 
However, once the physics of the stellar surface is mod- 
eled more realistically (e.g., without imposing an artifi- 
cial cooling layer to model radiation) it would be useful 
to produce synthetic surface maps and light curves that 
can be compared with observations. 

4. CONCLUSIONS 

The present work suggests that fully convective stars 
are capable of generating not only turbulent magnetic 
fields, but also strong large-scale fields that dominate 
the magnetic energy spectrum. In most of our models, 
the large-scale field has a strong quadrupolar component, 
in contrast to what is expected from mean-field theory 
for dynamo ac tion in thick shells and in full spheres 
l|Robertslll972|) . We have so far not seen evidence of 
magnetic cycles. The resolution of our models is still 
too low to be able to tell whether this type of magnetic 
field generation will continue to operate at much larger 
magnetic Reynolds numbers, but our results disprove the 
claim that a strong shear layer or a stably stratified core 
are necessary ingredients for the generation of large-scale 
magnetic fields. As one would expect in the absence of 
strong shear layers, the toroidal and poloidal components 
of the mean magnetic field are roughly comparable. 



Another important result concerns the self-consistently 
produced differential rotation. In our simulations, the 
angular velocity shows some tendency to be constant 
along cylinders, which is plausible for rapidly rotating 
stars. Whether or not this is realistic is difficult to say. 
Asteroseismology may in the future be able to reveal 
the internal angular velocity of stars, but at present the 
time coverage is still too short and incomplete. There 
is at least some hope of observing the surface differ- 
ential rotation, at least of sufficiently rapidly rotating 
fully convecti ve stars such as T Tauri sta rs, using sur- 
face imaging ijCollier Cameron et al.ll2004|) . This would 
be particularly interesting, given that our simulations 
predict a more slowly rotating equator. This behavior 
is opposite to that in the solar case. Thus far, the- 
ory i n terms of the A effect (e.g. iRudiger fc Hollerbachl 
|2004|) also tends to produce a faster equator, unless the 
turbulent motions possess a predominantly radial struc- 
ture (u rtimB 3> %j,rms)- In our case, however, there is 
strong meridional circulation, which, due to conservation 
of angular mome ntum, causes the outer l ayers to rotate 
more slowly; see Kitchatinov & Rudiecr (2004). Again, 
this result may no longer hold in real stars, because in 
our model the degree of stratification is far too low and 
the luminosity too high, so the convective velocities and 
meridional circulation tend to be exaggerated. 

Another reason for the slowly rotating equator could 
be connected with the outer boundary condition. In 
connection with geodynamo simulations there are indi- 
cations that a no-slip outer boundary condition (with 
respect to a rigidly rotating sp here) tends to produce a 
more slowly rotating equator ijChristensen et al1 fl999). 
In the limit of a short damping time, our effective outer 
boundary condition at r = R should indeed be closer 
to a no-slip condition than to a free-slip condition. On 
the other hand, in strongly magnetized stars the coro- 
nal magnetic field may enforce a rigidly rotating exte- 
rior, and hence produce conditions close to what is rep- 
resented by our model. 
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APPENDIX 
REFERENCE MODEL 

In order to specify the initial conditions and the gravity potential, we use a simple spherically symmetric, hydrostatic 
self-gravitating, isentropic model. The equations for this reference model are 

dm 



■ = 4irr g , 



dr 

dg Gm r g 2 ~^ 
dr jKr 2 

where m r denotes the total mass inside the sphere of radius r, together with the boundary conditions 

m r (0) = , g{0) = g c . 



(Al) 
(A2) 

(A3) 
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Here K = e 7S °/ Cp = const is the polytropic constant, relating pressure p and density g via p — Kg 1 , and sq is the 
constant value of entropy. The adiabatic exponent is 7 = 5/3, and thus our reference model is a polytropic model 
with polytro pic i ndex m _ 3/2. 

Equations (|A1|) and (|A2(I are integrated outwards, starting with certain values (g c , sq) for central density and entropy. 
As is common with polytropic models, the solution can have a surface (where g = p = T = 0) at some finite radius 
i? sur f, which must not be smaller than the desired stellar radius R. 

Varying the central values (g Cl sq), we can tune the reference model to match a given reference stellar radius R and 
total mass M. We choose the values R = 0.27 Rq = 1.9 x 10 8 m, and M = 0.21 M© = 4.2 x 10 29 kg, which correspond 
to an M5 dwarf with a luminosity of L = 0.008 L© = 3 x 10 24 W. 

Once the temperature and density profiles are known, one can calculate the approximate volume heating rate 
according to the formula 



H(g,T) « 4.9 x 10- [j^^^ (^j^^ . (A4) 

The T 53 dependence is an app roximation for the ppl chain of hydrogen burning near T c « 6 x 10 6 K; see Sec. 18.5.1 of 
iKippenhahn Sz Weigertl l)1990|) . We have used a Gaussian approximation to this (time-independent) radial dependence 
of H(r) for all simulations presented here, while adjusting the total luminosity by a multiplicative factor. Figure ^ 
shows the luminosity L r (r) as a function of radius. 

HIGH-ORDER UPWIND DERIVATIVES 

Convection simulations with high-order centered finite-difference schemes sometimes show a tendency to develop 
'wiggles' (Nyquist zigzag) in lag. This can be avoided by using a high-order upwind derivative operator, where the 
point furthest downstream is excluded from the stencil. We apply this technique only to the terms u • V In g and 
u • Vs. In the following we discuss the treatment in the x direction, but the treatment for the other directions is 
analogous. For u x > 0, we replace 

n , ~/-3 + 9/_ 2 - 45/_! + 45/i - 9/ 2 + h t , , fa 6 / (7) (&) ( 

-D C cnt,6 JO = = JO H ^40 ' 3 ^ 3 ' ^ ) 

by 

n , -2/_ 3 + 15/-2 - 6O/-1 + 20/q + 30/i - 3/ 2 ^ 5 / (6) fe) ^ t ^ 

Aip,5 JO = TTTTT = Jo ™ i X -3 < ?5 < %2 ■ (B2) 

60 ox 60 

Both formulae follow from Markoff's formula l|Abramowitz &: StegunllT984l §25.3.7). 

The difference between the sixth-order central and fifth-order upwind derivative is proportional to the sixth derivative 
operator 

n o , _ /-3 - 6/- 2 + 15/-i - 20/p + 15/x - 6/ 2 + h _ , (6 ) , Sx 2 f^(x ) 4 ^ 

-^ccnt^ Jo — — Jo ' | o'l^oa; ; , ^r>a; 

namely 

£>up,5 /o = £»ccnt, 6 /o - a 6x 5 D 6 Qcnt 2 , (B4) 

with a = 1/60. This allows us to represent the fifth-order upwind scheme in the advection term (for both signs of u x ) 
by sixth-order hyper-diffusion: 

- = -Uz/ccnt,6 + I ^ 5 /cont,2 ■ ( B5 ) 
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